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We demonstrate the synthesis of sparse sampling and machine learning to characterize and model 
complex, nonlinear dynamical systems over a range of bifurcation parameters. First, we construct 
modal libraries using the classical proper orthogonal decomposition to uncover dominant low-rank 
coherent structures. Here, nonlinear libraries are also constructed in order to take advantage of the 
discrete empirical interpolation method and projection that allows for the approximation of nonlinear 
terms in a low-dimensional way. The selected sampling points are shown to be nearly optimal 
sensing locations for characterizing the underlying dynamics, stability, and bifurcations of complex 
systems. The use of empirical interpolation points and sparse representation facilitate a family 
of local reduced-order models for each physical regime, rather than a higher-order global model, 
which has the benefit of physical interpretability of energy transfer between coherent structures. 

In particular, the discrete interpolation points and nonlinear modal libraries are used for sparse 
representation to classify the dynamic bifurcation regime in the complex Ginzburg-Landau equation. 

It is shown that nonlinear point measurements are more effective than linear measurements when 
sensor noise is present. 


I. INTRODUCTION 

The theoretical study of complex systems pervades 
the physical, biological and engineering sciences. Today, 
these studies are driven increasingly by computational 
simulations that are of growing complexity and dimen¬ 
sion due to numerical discretization schemes. Yet most 
dynamics of interest are known ultimately to be low¬ 
dimensional in nature (lj, thus contrasting, and in an¬ 
tithesis to, the high-dimensional nature of scientific com¬ 
puting. Reduced order models (ROMs) are of growing 
importance in scientific applications and computing as 
they help reduce the computational complexity and time 
needed to solve large-scale, complex systems [2|. Specif¬ 
ically, ROMs provide a principled approach to approx¬ 
imating high-dimensional spatio-temporal systems, typ¬ 
ically generated from numerical discretization, by low¬ 
dimensional subspaces that produce nearly identical in¬ 
put/output characteristics of the underlying nonlinear 
dynamical system. However, despite the significant re¬ 
duction in dimensionality, the complexity of evaluating 
higher-order nonlinear terms may remain as challenging 
as that of the original problem [3, 3 • The empirical in¬ 
terpolation method (EIM), and the simplified discrete 
empirical interpolation method (DEIM) for the proper 
orthogonal decomposition (POD) [g, 6|, overcome this 
difficulty by providing a computationally efficient method 
for discretely (sparsely) sampling and evaluating the non¬ 
linearity. These methods ensure that the computational 
complexity of ROMs scale favorably with the rank of the 
approximation, even with complex nonlinearities. 

An alternative computational strategy for handling the 
nonlinearity is based upon machine learning techniques 
whereby libraries of learned POD modes can be con¬ 
structed and inner products pre-computed for a num¬ 
ber of distinct dynamical regimes of the complex sys¬ 


tem [7I-IT0I]. This strategy also evokes the power of com¬ 
pressive sensing for efficiently identifying the active POD 
subspace necessary for a low-dimensional Galerkin-POD 
truncation Hi- In this manuscript, we combine the 
power of the DEIM with the library building strategy. 
Specifically, we show that building libraries that encode 
the nonlinearities allows one to (i) take advantage of 
DEIM to evaluate the nonlinearities, (ii) more robustly 
classify the dynamical regime the system is in, and (iii) 
identify the discrete and optimal sensor locations to eval¬ 
uate a nonlinear model reduction. We demonstrate the 
full integration of the methods on a canonical model of 
mathematical physics and nonlinear science, the cubic- 
quintic Ginzburg-Landau (CQGLE) equation. 


A. Dimensionality Reduction 

Although a variety of dimensionality-reduction tech¬ 
niques exist, the ROM methodology considered here is 
based upon the proper orthogonal decomposition Q. 
The POD method is ubiquitous in the dimensionality 
reduction of physical systems. It is alternatively re¬ 
ferred to as principal components analysis (PCA) EL 
the Karhunen Loeve (KL) decomposition, empirical or¬ 
thogonal functions (EOF) EL or the Hotelling trans¬ 
form IH, 0 . Snapshots (measurements) of many 
complex system often exhibit low-dimensional phenom¬ 
ena (l|, so that the majority of variance/energy is con¬ 
tained in a few modes computed from a singular value 
decomposition (SVD). For such a case, the POD basis 
is typically truncated at a pre-determined cut-off value, 
such as when the modal basis contain 99% of the vari¬ 
ance, so that only the first r-modes (r-rank truncation) 
are kept. There are numerous additional criteria for 
the truncation cut-off, and recent results derive a hard- 
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threshold value for truncation that is optimal for sys¬ 
tems with well-characterized noise [15]. The SVD acts 
as a filter, and so often the truncated modes correspond 
to random fluctuations and disturbances. If the data 
considered is generated by a dynamical system (nonlin¬ 
ear system of ordinary differential equations of order n), 
it is then possible to substitute the truncated POD ex¬ 
pansion into the governing equation and obtain Galerkin 
projected dynamics on the rank-r basis modes @, 113 . 
Recall that we are assuming that the complex systems 
under consideration exhibit low-dimensional attractors, 
thus the Galerkin truncation with only a few modes 
should provide an accurate prediction of the evolution 
of the system. Note that it has also been shown re¬ 
cently that it is possible to obtain a sketched -SVD by 
randomly pro jecting the data initially and then comput¬ 
ing the SVD [lMl. 


B. Sparse Sampling 

EIM has been developed for the purpose of efficiently 
managing the computation of the nonlinearity in dimen¬ 
sionality reduction schemes, with DEIM specifically tai¬ 
lored to POD with Galerkin projection. Indeed, DEIM 
approximates the nonlinearity by using a small, discrete 
sampling of points that are determined in an algorithmic 
way. This ensures that the computational cost of evalu¬ 
ating the nonlinearity remains proportional to the rank 
of the reduced POD basis. As an example, consider the 
case of an r-mode POD-Galerkin truncation. A simple 
cubic nonlinearity requires that the POD-Galerkin ap¬ 
proximation be cubed, resulting in r 3 operations to eval¬ 
uate the nonlinear term. DEIM approximates the cubic 
nonlinearity by using 0{r) discrete sample points of the 
nonlinearity, thus preserving a low-dimensional (O(r)) 
computation, as desired. The DEIM approach combines 
projection with interpolation. Specifically, DEIM uses 
selected interpolation indices to specify an interpolation- 
based projection for a nearly optimal £2 subspace approx¬ 
imating the nonlinearity. EIM/DEIM are not the only 
methods developed to reduce the complexity of evaluat¬ 
ing nonlinear terms, see for instance the missing point 
estimation (MPE) [lj| or gappy POD flil [22| methods. 
However, they have been successful in a large number 
of diverse applications and models [J. In any case, the 
MPE, gappy POD, and EIM/DEIM use a small selected 
set of spatial grid points to avoid evaluation of the expen¬ 
sive inner products required to evaluate nonlinear terms. 

The discrete sampling points given by DEIM to evalu¬ 
ate the nonlinearity get a new interpretation in the cur¬ 
rent work. Specifically, we show them to be the nearly 
optimal locations for placing sensors in the complex sys¬ 
tem in order to (i) determine the dynamic regime of the 
system, (ii) reconstruct the current state of the system, 
and (iii) produce a POD-Galerkin prediction (nonlinear 
model reduction) of the future state of the system. Such 
tasks are accomplished by using ideas of sparse repre¬ 


sentation [23} and compressive sensing [24M3i|. In par¬ 
ticular, the theory of compressive sensing shows that a 
small number of measurements are sufficient to perform 
a reconstruction provided there exists a sparse represen¬ 
tation (or basis) of the data. Sparsity techniques have 
also been shown to be highly effective for numerical so¬ 
lution schemes (H, [33j. In our case, the sparse basis is 
generated from a library learning procedure. More than 
that, however, we also build libraries of the nonlineari¬ 
ties , thus pre-computing the low-dimensional structures 
observed in the different dynamical states of the complex 
system. This allows for more robust dynamical classifica¬ 
tion as well as allowing easy evaluation of the nonlinear 
terms through DEIM. The combination of library build¬ 
ing, compressive sensing and DEIM is demonstrated to 
be a highly effective and intuitively appealing methodol¬ 
ogy for scientific computing applications. It further high¬ 
lights the need in modern scientific computing of complex 
systems to integrate a variety of data-driven modeling 
strategies, many of which are being developed under the 
aegis of machine learning, in order to most efficiently 
simulate large-scale systems. 


C. Physical Interpretation 

The ideas presented here are more than just numerical 
efficiencies. Indeed, the methodology identifies the un¬ 
derlying modal structures that drive the dynamics of the 
complex system, thus helping to understand the funda¬ 
mental interactions and physics of the system. Through¬ 
out the development of 20th-century physics and engi¬ 
neering sciences, the understanding of many canonical 
problems has been driven by recasting the problem into 
its natural basis (mode) set. The majority of classical 
problems from mathematical physics are linear Sturm- 
Liouville problems whose ideal modal representations are 
generated from eigenfunction decompositions, i.e. special 
functions. In quantum mechanics, for instance, Gauss- 
Hermite (denoted by H n {x)) polynomials are the natural 
basis elements for understanding the harmonic oscilla¬ 
tor. Likewise, spherical harmonics (denoted by Y[ m (0, <p)) 
are critical in the computation of atomic orbital elec¬ 
tron configurations as well as in representation of grav¬ 
itational fields, the magnetic fields of planetary bodies 
and stars, and characterization of the cosmic microwave 
background radiation. 

For modern complex systems, nonlinearity plays a 
dominant role and shapes the underlying modes, thus 
necessitating a new approach, such as that presented 
here, for extracting these critical spatio-temporal struc¬ 
tures. Remarkably, although nonlinearity creates new 
modal structures, it does not destroy the underlying low¬ 
dimensional nature of the dynamics. Distinct physical 
regimes may be obtained by varying bifurcation parame¬ 
ters, and these regimes will typically have different local 
bases and physical interactions. Instead of developing a 
global interpolated model, which may obscure these dis- 
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tinct physical mechanisms, we advocate a hierarchy of 
models along with sparse sampling and machine learning 
to classify and characterize the system parameters from a 
few online measurements. Methods that take advantage 
of such underlying structure are critical for developing 
theoretical understanding and garnering insight into the 
fundamental interactions of the physical, engineering and 
biological systems under consideration. 

The paper is outlined as follows. In Sec. El an overview 
of the mathematical framework of the POD method and 
the DEIM is given. This is followed up in Sec. IIIII with 
an introduction of the nonlinear dynamical system, i.e. 
the cubic-quintic Ginzburg-Landau equation, where the 
methods proposed here will be applied. The library 
building procedure that encodes the various dynamical 
regimes of our model equation are discussed in Sec. IIVI 
Once the libraries are constructed, DEIM points, or sen¬ 
sor locations, are computed in Sec. IVland their ability to 
classify dynamical regimes is evaluated in Sec. [VI] The 
reconstruction of the dynamics and future state projec¬ 
tion is discussed in Sec. I VIII A summary of our findings 
and an outlook on the method is given in the concluding 
Sec. rvTTTl 

II. BACKGROUND FOR MODEL REDUCTION 

Our innovations are built upon two key methods which 
are used for model reduction and approximating nonlin¬ 
ear dynamical systems. The first approach is the well- 
known POD-Galerkin method, which is used to reduce 
the dimension of systems in a principled way. However, 
computing the form of the nonlinearity in the reduced- 
order system is an expensive offline computation, as in¬ 
ner products of the full high-dimensional system must 
still be computed. Online evaluation of the nonlinear 
terms in the reduced order model may remain expensive, 
as these typically involve dense matrix or tensor opera¬ 
tions of the same order as the degree of nonlinearity. The 
second approach highlighted is the DEIM algorithm 
which reduces the complexity of evaluating the nonlinear 
terms. In particular, it gives a principled way to sparsely 
sample the nonlinearity in order to approximate the non¬ 
linear terms in a low-dimensional way. 

A. POD 

Consider a high-dimensional system of nonlinear differ¬ 
ential equations that can arise, for example, from the fi¬ 
nite difference discretization of a partial differential equa¬ 
tion: 

^ = Lu(t) + N(u(t)), (1) 

where u(t) = [ui(t) u 2 (t) ••• u n (t)] T £ R" and n 3> 
1. Typically under discretization of a single spa¬ 
tial variable, Uj{t) = u(xj : t) is the value of the 


field of interest at the spatial location Xj. The lin¬ 
ear part of the dynamics is given by I £ R raxrl 
and the nonlinear terms are in the vector N(u(t)) = 
[7Vi(u(t)) N 2 (u(t)) ■■■ N n (u(t))] T £ R". The non¬ 

linear function is evaluated component-wise at the n spa¬ 
tial grid points used for discretization. 

For achieving high accuracy solutions, n is typi¬ 
cally required to be a very large number, thus mak¬ 
ing the computation of the solution expensive and/or 
intractable. The POD-Galerkin method is a princi¬ 
pled dimensionality-reduction scheme that approximates 
the function u(i) with rank-r optimal basis functions 
where r <C n. These optimal basis functions are com¬ 
puted from a singular value decomposition of a series 
of temporal snapshots of the complex system. Specif¬ 
ically, suppose snapshots of the state, u (tj) with j = 
1 , 2 , ,p, are collected. The snapshot matrix X = 

[u(fi) u(t 2 ) ••• u(t p )] £ R" xp is constructed and the 
SVD of X is computed: X = $£W*. The r-dimensional 
basis for optimally approximating u(£) is given by the 
first r columns of matrix 3?, denoted by 4> r . Thus the 
POD-Galerkin approximation is given by 

u(i) w $ r a(i) (2) 

where a(t) £ R r is the time-dependent coefficient vector 
and r <C n. Plugging this modal expansion into the gov¬ 
erning equation © and applying orthogonality (multi¬ 
plying by <!>(/) gives the dimensionally reduced evolution 

^ = *?L* r a.(t) + *JTlV(* r a(i)). (3) 

By solving this system of much smaller dimension, the 
solution of a high-dimensional complex system can be 
approximated. 

This standard POD procedure [gj has been a ubiqui¬ 
tous algorithm in the reduced order modeling commu¬ 
nity. However, it also helps illustrate the need for inno¬ 
vations such as DEIM, Gappy POD and/or MPE. Con¬ 
sider the nonlinear component of the low-dimensional 
evolution ©: <frj r ./V(<l> r a(f)). For a simple nonlinear¬ 
ity such as N(u(x,t)) = u(x,t) 3 , consider its impact 
on a spatially-discretized, two-mode POD expansion: 
u(x : t) = a\{t)4>\{x)+a 2 {t)(j> 2 {x). The algorithm for com¬ 
puting the nonlinearity would require the evaluation: 

u(x, t ) 3 = a\(j)\ + 3a^02</i</2 + • (4) 

The dynamics of ai(t) and 02 (f) would then be computed 
by projecting onto the low-dimensional basis set by tak¬ 
ing the inner product of this nonlinear term with respect 
to both (fii and <j) 2 . Thus the number of computations not 
only doubles, but the inner products must be computed 
with the n-dimensional vectors. Methods such as DEIM 
overcome this high-dimensional computation and instead 
produce an 0{r) dimensional evaluation of the nonlinear 
terms. 
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TABLE I: DEIM algorithm for finding approximation basis for the nonlinearity and its interpolation indices. 


DEIM algorithm 

Basis 

• collect data, construct snapshot matrix 

X = [u(ti) u(t 2 ) ■ • ■ u(f p )] 

• construct nonlinear snapshot matrix 

N = [JV(u(i 1 )) JV(u(i 2 )) ••• N{u{t p ))\ 

• singular value decomposition of N 

N = HEjvW^r 

• construct approximating basis (first m columns) 

= [£l ^2 ’ ’ ’ £m] 

Interpolation Indices (Iteration Loop) 

• choose the first index (initialization) 

[p, 7i] = max |£i | 

• approximate by £1,..., £7-1 at indices 71, ...,7^-1 

Solve for c: P J = P J 3 j ic with P = [e 71 ■ • • e 7 , 1 ] 

• select 7j and loop ( j = 2 , 3,..., m) 

[p,7j] = max |£j - Hj-ic| 


B. DEIM 

As outlined in the previous section, the shortcomings 
of the POD method are generally due to the evaluation of 
the nonlinear term N(<& r a.(t)). To avoid this difficulty, 
the DEIM approximates N = iV(3> r a(t)) through pro¬ 
jection and interpolation instead of evaluating it directly. 
Specifically, a low-rank representation of the nonlinearity 
is computed from the singular value decomposition 

N = EEatW^ (5) 

where the matrix 3 contains the optimal (in an £2 sense) 
basis set for spanning the nonlinearity. Specifically, we 
consider the rank-m basis set 3 m = [£1 £2 • • • £ m ] that 
approximates the nonlinear function (m <C n and m ~ r). 
The approximation to the nonlinearity N is given by: 

N « S ro c (t) (6) 

where c(t) is similar to a (t) in ©. Since this is a highly 
overdetermined system, a suitable vector c (t) can be 
found by selecting only m rows of the system. The DEIM 
algorithm was specifically developed to identify which m 
rows to evaluate. 

The DEIM algorithm begins by considering the vec¬ 
tors e 7j £ R” which are the 77 -th column of the n di¬ 
mensional identity matrix. We can then construct the 
projection matrix P = [e 7l e 72 • • • e 7m ] which is cho¬ 
sen so that P 2 3 m is nonsingular. Then c(t) is uniquely 
defined from P T N = P 2 3 m c(t), and thus, 

N « H m (P T H ro )- 1 P T N. (7) 

The tremendous advantage of this result for nonlinear 
model reduction is that the term P T N requires evalua¬ 
tion of nonlinearity only at m indices, where m <C n. The 
DEIM further proposes a principled method for choosing 
the basis vectors and indices 7 j. The DEIM algorithm, 
which is based upon a greedy-like search, is detailed in Q 
and further demonstrated in Table [0 


C. Application to ROMs 

POD and DEIM provide a number of advantages for 
nonlinear model reduction of complex systems. POD pro¬ 
vides a principled way to construct an r-dimensional sub¬ 
space 3> r characterizing the dynamics. DEIM augments 
the POD method by providing a method to evaluate the 
problematic nonlinear terms using an m-dimensional sub¬ 
space S m that represents the nonlinearity. Thus a small 
number of points, specifically m, can be sampled to ap¬ 
proximate the nonlinear terms in the ROM. 

The method proposed here capitalizes on these meth¬ 
ods by building low-dimensional libraries associated with 
the full complex system dynamics as well as the specific 
nonlinearities. Moreover, the sparse measurement loca¬ 
tions computed by DEIM are found to be nearly optimal 
for sensor placement. Such sensors, as will be shown in 
what follows, can be used with sparse representation and 
compressive sensing to (i) identify dynamical regimes, (ii) 
reconstruct the full state of the system, and (iii) provide 
an efficient nonlinear model reduction and POD-Galerkin 
prediction for the future state. Moreover, we show that 
nonlinear measurements of the dynamical system can be 
much more robust to noise for accomplishing the above 
tasks. 

III. MODEL PROBLEM 

One of the canonical nonlinear PDEs in mathematical 
physics and pattern forming systems is the Ginzburg- 
Landau (GL) equation and its many-variants jl[. It has 
been used to model a variety of physical systems from 
condensed matter to biological waves. Here we consider 
a variant of the GL equation arising in mode-locked laser 
theory that has cubic and quintic nonlinear terms and a 
fourth-order derivative (34j : 

iUt -j- iT ^ U X x ikU x xxx + (1 i/i)|t/j U 

+{v-ie)\U\ 4 U -^11 = 0 , (8) 

where U(x,t) is a complex valued function of space and 
time. Under discretization of the spatial variable, U ( x , t) 
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T 

K 

d 

V 

e 

7 

description 

Pi 

-0.3 

-0.05 

1.45 

0 

-0.1 

-0.5 

3-hump, localized 

02 

-0.3 

-0.05 

1.4 

0 

-0.1 

-0.5 

localized, side lobes 

p3 

0.08 

0 

0.66 

-0.1 

-0.1 

-0.1 

breather 

Pi 

0.125 

0 

1 

-0.6 

-0.1 

-0.1 

exploding soliton 

05 

0.08 

-0.05 

0.6 

-0.1 

-0.1 

-0.1 

fat soliton 

P6 

0.08 

-0.05 

0.5 

-0.1 

-0.1 

-0.1 

dissipative soliton 


TABLE II: Values of the parameters from equation (j8j that 
lead to six distinct dynamical regimes. To exemplify our al¬ 
gorithm, the first, third and fifth regimes will be discussed in 
this paper. 


becomes a vector u with n components, i.e. u j(t) = 
U ( Xj,t ) with j = 1 , 2 , • • • n. 

An efficient and exponentially accurate numerical so¬ 
lution to (0 can be found using standard spectral 
methods |l£fl. Specifically, the equation is solved by 
Fourier transforming in the spatial dimension and then 
time-stepping with an adaptive 4th-order Runge-Kutta 
method. The extent of the spatial domain is x £ [—20, 20] 
with n = 1024 discretized points. Note that in what 
follows, the indices for evaluation of the nonlinear term 
correspond to the collocation points away from the cen¬ 
ter spatial point of the computational domain 2513 = 0 . 
Here, we allow the parameters /3 = (r, k, h, v, e, 7 ) to vary 
in order to discover various dynamical regimes that ex¬ 
hibit low-rank structure and stable attractors. Table El 
shows six different parameter regimes that have unique 
low-dimensional attractors (see |9[). The evolution of the 
system for parameter regimes p \, P 3 and /3s is illustrated 
in Fig. [Q Such stereotypical low-dimensional behaviors, 
which are commonly observed in pattern forming sys¬ 
tems [li], will serve as the basis for our library building 
methodology, especially in regards to using a small num¬ 
ber of measurements to identify the /3j regime, recon¬ 
struct the solution, and project a future state. Although 
our results are demonstrated on this specific PDE, the 
methodology is quite general. 


IV. LIBRARIES 

As can be seen from Fig. [T| and Table El generic initial 
conditions evolve towards a variety of low-dimensional 
attractors. This suggests that each dynamic regime, with 
a given f3j, can be approximated by a small number of 
modes via a POD reduction. These modes will constitute 
our library modes in what follows. For each of the six 
regimes /3j in Table El we build a library of POD modes. 
The number of POD modes r is selected to capture 99% 
of the total variance (energy). For the /3i, /3 2 , /3s and /3q 
regimes, only a single mode is required so that r = 1. For 
the P 3 regime r = 6 , whereas for the 04 regime, r = 14 
in order to capture the fluctuations observed. Figure 
Ha) illustrates the library POD modes in differing colors 
for all of the y 3j regimes except /34. The exclusion of 



FIG. 1: Evolution dynamics of 0 for the parameter regimes 
/3i, /?3 and /3s over the time interval t £ [0,40]. The initial 
transients are quickly attenuated away, leaving the stable at¬ 
tractor for the given /3j regime. Sampling of the dynamics for 
library building occurs once the transients have decayed. 


the P 4 modes in this visualization is simply due to the 
large number (r = 14) necessary in comparison to the 
other dynamical regimes. As illustrated in Fig.[ 2 j library 
building is the first step in a training module aimed at 
learning the low-rank dynamical behavior of a complex 
system. 

In practice, a dynamical system such as 0 may change 
over time due to evolution or modulation of the param¬ 
eters /3j. Thus the dynamics may evolve from one at¬ 
tractor to another with some prescribed transition time 
(typically on the order of 0(1) time for 0). One of 
the primary goals of this and previous [35j work is to 
find optimal and sparse sensor locations whereby limited 
measurements of the system are taken in order to clas¬ 
sify the dynamical regime. Interestingly, the previous 
efforts [j] used expert-in-the-loop knowledge to help se¬ 
lect the optimal measurement positions. For the simple 
model considered here, such expert knowledge can be ac¬ 
quired from familiarity with the POD library modes and 
considering locations of maximal variance. However, for 
a more general system, this is a difficult task that could 
greatly benefit from a more principled mathematical ap¬ 
proach. The DEIM algorithm will provide this approach. 
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1. Data Collection 2. POD Modes + 
_ Nonlinear POD 


I. Training Module 


Complex System 
u = Lu + N(u) 


II. Execution Module 




3. DEIM Algorithm 
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1. Sparse 

2. Sparse 3. 

Reconstruction** 4. Galerkin-POD** 

Measurement 

Classification 

( /3j from step 2) Projection 


c = argmin c / ||c'||i 


u = Pu 

subject < 
to* 

u — 

N 3 = P\P 3 C U = ^L,/3 j (P^L,/3 J ) t U U=$ I/ ^.a(t) 

N 5 = P^sc 

, Nm = P’FjviC 


*Note: one of the above 

** Note: &l,i 3 » contained in 'b l 


FIG. 2: Training and execution modules for the library learning and sensor location optimization with DEIM. The training 
module samples the various dynamical regimes (/3i, /?2 , ■ • • ,/3j) through snapshots. For each dynamical regime, low-rank 
libraries are constructed for the nonlinearities of the complex system (<&£«., .. <&5 , &NL,Pj)- The DEIM algorithm is 

then used to select sparse sampling locations and construct the projection matrix P. The execution module uses the sampling 
locations to classify the dynamical regime /3j of the complex system, reconstruct its full state (u = 3?L,/3j (PH>L,£j)'u), and 
provide a low-rank Galerkin-POD approximation for its future (u = $g^a(t)). Note that (P^l,/ 9 j)* denotes the Moore- 
Penrose pseudo-inverse of (P4?L,/3j)- 


Moreover, as required by DEIM, we also build low-rank 
libraries for the cubic and quintic terms associated with 
the dynamical regimes /3j. In doing so, we not only find 
nearly optimal sensor locations, but we also circumvent 
the computational difficulties of the POD in evaluating 
the nonlinear terms. 

To library build, consider the following linear and non¬ 
linear functions associated with the governing equations 
© for a given parameter regime (3j: 

N l (U) = U (9a) 

N 3 {U) = \U\ 2 U (9b) 

N 5 {U) = \U\ A U (9c) 

N nl {U) =(i + p)\U\ 2 U + 0 iv + e)\U\ A U , (9d) 

where the second and third terms are the standard cu¬ 
bic and quintic nonlinearities of (JSJl and the last term 


enforces their prescribed relative weighting. 

Associated with each nonlinearity <[9|) are a set of mea¬ 
surements and snapshot matrices. For a snapshot matrix 
sampled at p temporal locations [u 3 u 2 • • ■ u p ] £ R raxp , 
we can construct the nonlinear R™ xp snapshot matrices 

N l = [ui u 2 • • • Up] (10a) 

N 3 = [7V 3 (ui) N 3 (u 2 ) ••• Ar 3 (u p )] (10b) 

N 5 = [iV 5 (ui) IV5(u 2 ) ••• JV 5 (up)] (10c) 

Natl = [N nl ( u 3 ) N nl ( u 2 ) ••• N NL (u p )]. (lOd) 

The singular value decomposition of these matrices will 
give a basis for approximation of each of the nonlineari¬ 
ties for a given /3j as well as the standard snapshot matrix 
of POD. Specifically, the SVD gives the library of modes: 

3 * 3 , Pj, ^ 5 , 0 j and ^NL,Pj (See Fig. [2]). 

The POD modes can be arranged in a collection of 
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modes of the full system 



nonlinear cubic modes 



nonlinear quintic modes 



FIG. 3: Library modes for (a) the full system, (b) the cubic 
nonlinearity, and (c) the quintic nonlinearity. The modes are 
color coded by their dynamical regime from /ft to /ft as given 
in Table [TT] The rank-r for each library is chosen by selecting 
the modes that comprise 99% of the total variance for a given 
dynamical regime. 


library elements, \& 3 , IE's or \HjVL, by concatenating 
the POD modes from each of the different Bj regimes. 
Thus the construction of multiple libraries would take 
the form 


IE'l = [$L,ft ®L,f 3 2 ■ 


(11a) 

V E r 3 = [3ft,ft ‘ ’ 

’ ^3, ft] 

(lib) 

^5 = [$5,ft $5, ft ‘ ' 

• $5,/3 6 ] 

(11c) 

&NL = l&NLj 3i &NL,/3 2 ' ' ' &NL,/3 e \ ■ 

(lid) 


The number of basis elements (rank) for the cubic and 
quintic terms in a given POD library coincides with the 
rank r required for each /3j, i.e. r = m. Note that the li¬ 
brary ’S'l is the library containing the POD modes used 
for POD-Galerkin projections of the future state. It is 
also the only library constructed in previous work [3, Q. 
Figure [3Kb, c) shows the cubic and quintic library modes 
for ([8|. They can be compared to the standard POD 
modes illustrated in Fig. [dja). Although the modes look 
quite similar, we will show that the classification can be 
improved with the nonlinear libraries. Further, evalua¬ 
tion of the nonlinearities through DEIM now remains a 
low-order computation. 


V. DEIM FOR SENSOR LOCATIONS 

The idea of using a limited (sparse) number of sensors 
to characterize the dynamics has previously been con¬ 
sidered in However, no algorithm was specified 


to determine the best locations for the sensors, although 
optimal sensor placement has been investigated in the 
context of categorical decisions (35j. Indeed, the pre¬ 
vious work relied on expert-in-the-loop selection of the 
sensors in order to classify the dynamics. Interestingly, 
the DEIM algorithm gives a principled way to discretely 
and sparsely sample the nonlinearity in order to evaluate 
the various inner products for a POD reduction. This 
begs the question: would these same DEIM spatial sam¬ 
pling locations make good sensor locations for classifica¬ 
tion and reconstruction? Since the interpolation indices 
from the DEIM algorithm Q correspond to the entries 
with largest magnitude of the residual error between the 
chosen basis and its approximation at each step (see last 
line of the table H|, it becomes interesting to see what the 
classification/reconstruction will be if we pick these loca¬ 
tions for sensors. As demonstrated in Fig. [2j determining 
the sensor locations is part of a training module. 

We apply the DEIM algorithm outlined in Table U 
on the nonlinear POD (SVD) library modes (HG, ^5 
or nl ) computed from m and m ■ The applica¬ 
tion of the algorithm yields DEIM interpolation loca¬ 
tions which we will call our sensor locations. Note that 
the indices indicate the distance away from the center 
of the computational grid. Thus xq = 0, x±\ = dx, 
x ±2 = 2dx , etc. Or more generally, the index n corre¬ 
sponds to x n = ndx. Thus the indices depend on the 
specific discretization of the domain. Sensor locations 
are computed for each of the nonlinearities: < I > 5 ,ft 

and 3>ArL,ft for j = 1,2,3. Each dynamical regime /3j 
and nonlinear library gives a unique set of sensor loca¬ 
tions. Our goal is to evaluate the placement of 3 sensors. 
Table IIIII and its accompanying figure gives a vector of 
the indices for the locations x^. of the 3 sensors found 
for three regimes /ft, /ft and /3s using the libraries ^ 3 ,ft, 
‘E’s,ft and ‘E’Afftft. Also represented are the 3 sensor lo¬ 
cations when all three /ft regimes are combined into a sin¬ 
gle library, i.e. the best sensor locations for the combined 
dynamic library is identified. This regime is represented 
in Table Hill by x^. 

Application of the DEIM algorithm results in the mea¬ 
surement matrix P of ©• For 3 sensors, generically it 
takes the form 


P = 


1 0 ••• ••• 0 

0 ••• 0 1 0 . 0 

0 ••• •••010 ■•• 0 


( 12 ) 


where the specific columns containing the nonzero entries 
are given by the indices found from DEIM and shown in 
Table Hm More precisely, this matrix is exactly the output 
of the DEIM algorithm. In our scenario, the construction 
of the P matrix is made for each nonlinearity as well 
as for each dynamical regime Bj . This gives the nearly 
optimal sensor locations for the sparse sensing scheme 
presented in the next section. Figure [4] illustrates the 
locations of the sensors and the value of library modes at 
the prescribed locations for both the cubic and quintic 
nonlinearities. 

















































DEIM interpolation indices 

Pi p3 Ps ball ••• $L 




FIG. 4: Location of indices determined by DEIM for the nonlinear libraries \U\ 3 , |t/| 5 and N(U). The spatial domain x £ 
[—20, 20] is discretized on a periodic domain with n = 1024 points. The center point of the domain corresponds to a;(0) = 0. 
The index values are the number of grid points ndx away from the center grid point, e.g. *(5) = 5 dx. The left grid shows 
the location of the DEIM indices (black boxes) determined by the algorithm in Table [T] for the regimes Pi, Pz and /3g as well 
as the combination of all three regimes together p a ii- The middle panel shows the library mode &l,Pi (laid out vertically) 
as a function of the spatial variable x(n). Indicated on this transverse mode are the measurement locations for the different 
DEIM nonlinearities and Pj regimes. The right two panels show the Pi, P 3 and Ps modes with the black lines indicating the 
measurement locations for n = 0, 6 and 13. This allows one to visualize where the measurement occur on the mode structures. 



Cubic 

\U\ 2 U 

Quintic 

\U\ 4 U 

Nonlinear 

N(U) 

Sensor 

XSl 

X <S 3 

x /3 5 


XSi 

X <S 3 

x /3 5 

X /3 all 

x /3l 

x /3 3 

x /3 5 

x /3„.n 

one 

0 

0 

0 

0 

0 

0 

0 

0 

0 

9 

0 

0 

two 

5 

15 

12 

6 

4 

13 

10 

6 

6 

21 

6 

6 

three 

13 

26 

17 

22 

13 

23 

15 

20 

13 

32 

15 

13 


TABLE III: Summary of sensor location vectors (indices for 
evaluation) from the DEIM algorithm. The table summarizes 
the findings from Fig. [4] giving precise grid cells to be used in 
evaluating the nonlinear inner products in the Galerkin-POD 
approximation. 


VI. CLASSIFICATION 

Our goal is to make use of recent innovations in sparse 
sampling and compressive sensing [2~fl |7TTI| for character¬ 
izing the complex system 0 - 0 - Specifically, we wish to 
use a limited number of sensors for classifying the dy¬ 
namical regime of the system. With this classification, 
a reconstruction of the full state space can be accom¬ 
plished and a POD-Galerkin prediction can be computed 
for its future. In general, if we have a sparse measure¬ 
ment u £ R 9 , where q is the number of measurements, 
then 

u = Pu, (13) 


where u is the full state vector and P is the sampling ma¬ 
trix determined by DEIM given by ( 1121 ) . In the previous 
section, we constructed the matrix P for q = 3. 

The full state vector u can be approximated with the 
POD library modes (u = \Flc), therefore 

u = P’S'lc, (14) 

where tF l is the low-rank matrix whose columns are POD 
basis vectors concatenated across all f) regimes and c 
is the coefficient vector giving the projection of u onto 
these POD modes. If P\Fl obeys the restricted isometry 
property (36] and u is sufficiently sparse in iFl, then it 
is possible to solve the highly-underdetermined system 
m with the sparsest vector c. Mathematically, this is 
equivalent to the optimization problem 

c = min ||c'||o, subject to u = P’J'lc. 

c' 

Minimizing the Iq norm is computationally an np-hard 
problem. However, It has been proven that under cer¬ 
tain conditions, a sparse solution of equation ( 1141 ) can be 
found by minimizing the 1 1 norm instead [25l. |27| so that 

c = argmin ||c'||i, subject to u = P\&£C. ( 15 ) 

c' 

The last equation can be solved through standard con¬ 
vex optimization methods such as the CVX package for 
Matlab. 
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cubic 



mode number n 


quintic 



0 5 10 15 20 25 

mode number n 


FIG. 5: The values of the 24 x 1 projection vector c from 
solving using a cubic measurement u 3 = |u| 2 u and the cu¬ 
bic library 'I'3 in (I16h l. The three panels show the dominant 
vector component to be in the /3i, ft 3 and /3s regime respec¬ 
tively, thus showing that it correctly identifies each dynami¬ 
cal regime from 3 measurement locations. The values of the 
colored circles correspond to the expression strength of the 
different library elements of Fig. [3] 


FIG. 6: The values of the 24 x 1 projection vector c from solv¬ 
ing using a quintic measurement U5 = |u| 4 u and the quintic 
library It's in (116b l. The three panels show the dominant 
vector component to be in the /3i, /J3 and /3s regime respec¬ 
tively, thus showing again that nonlinear measurements cor¬ 
rectly identify each dynamical regime from 3 measurement 
locations. The values of the colored circles correspond to the 
expression strength of the different library elements of Fig. [3] 


To classify the dynamical regime from limited measure¬ 
ments u (specifically 3 spatial measurements), we use the 
sensor locations matrix P found from DEIM on the non¬ 
linear libraries. Here, the sensor locations used for P are 
from all the library elements combined and the nonlin¬ 
earity N(U) (See the last column in Tabic IIIII remarked 
with red boxes), i.e. n = 0,6 and 13. Suppose we have a 
linear measurement u, then we can construct the vectors 
u 3 = |u| 2 u and U 5 = |u| 4 u and classify them using the 
nonlinear libraries. Specifically, the nonlinear classifica¬ 
tion is accomplished with: 

c 3 = argmin 11C 3 11 1 , subject to u 3 = P\F 3 c 3 (16a) 

C 3 

C 5 = argmin ||cg||i, subject to u 5 = P\F, 5 C 5 . (16b) 

c 5 

Figures [5] and [HI show the coefficient vectors c 3 and c 5 
respectively for measurements performed in the /%, /? 3 
and /3s regimes. The vectors c 3 and c 5 clearly act as 
accurate indicator functions for the dynamical regime. 
Indeed, the DEIM algorithm for sensor location does as 
well as expert-in-the-loop selections , but requires no 
extensive and pre-existing knowledge about the dynam¬ 
ical libraries. We can also make a categorical decision, 
with similar results, about the dynamical regime the dy¬ 
namics belongs to by computing error of projection onto 
a given library and considering which has the smallest 
error. This is the same as sparse representation used for 
image classification [23j. 

The above analysis assumes that there is no noise in 
the measurements or the system itself. However, most 


sensors are subject to noise fluctuations which can impact 
the ability of a scheme such as this to correctly identify 
fjj. As a consequence, we also perform the classification 
task with noisy data. First, assume that we collect linear 
measurements which have additive noise. Denote this 
data by 

u = u + A/”(0 ,ct 2 ) (17) 

where Af(Q, cr 2 ) is a Gaussian distributed noise term with 
variance a 1 . 

In order to evaluate the classification, we need to once 
again compute the nonlinear terms and run the optimiza¬ 
tion algorithm for computing the library coefficients and 
the associated dynamical regime. The statistical result 
for 400 trials when a = 0.2 is shown in Table El One 
can see that the noise introduces misclassification errors 
to the original 100% accurate classification scheme. How¬ 
ever, multiple measurements still give an accurate clas¬ 
sification overall with the exception of using the quintic 
library in the /? 3 regime. 

Interestingly, if nonlinear measurements are consid¬ 
ered, then the results can improve drastically. For in¬ 
stance, in optics, measurements are made of the intensity 
of the field rather than the field itself. This represents a 
simple form of a nonlinear measurement. Thus consider 
the nonlinear measurements subject to noise: 

u 3 = |u| 2 u + A/"(0,cr 2 ) (18a) 

u 5 = |u| 4 u + Af(Q, cr 2 ). (18b) 

The classification results for this case are also shown in 
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0i regime 

01 

$2 

03 

04 

05 

06 

|u| ‘‘u 

98.75 

0 

1.25 

0 

0 

0 

|m| 4 m 

91 

6.5 

2.5 

0 

0 

0 

U3 

100 

0 

0 

0 

0 

0 

u 5 

100 

0 

0 

0 

0 

0 








03 regime 

01 

P‘2 

03 

04 

Pb 

06 

\u\ z u 

2.5 

0 

61.75 

18 

17.5 

0.25 

|m| 4 m 

5.5 

0 

38 

34.5 

21.75 

0.25 

U3 

0 

0 

100 

0 

0 

0 

«5 

0 

0 

100 

0 

0 

0 








0 5 regime 

0i 

P2 

03 

04 

P5 

06 

\u\‘‘u 

5.25 

0.75 

7.5 

5 

62 

19.5 

|m| 4 m 

6.75 

2 

6.25 

2.5 

61.25 

21.25 

«3 

0 

0 

0 

0 

100 

0 

«5 

0 

0 

0 

0 

100 

0 


TABLE IV: Classification accuracy with noisy measurements 
(<r = 0.2) using 400 realizations in the 0\, 03 and 05 regimes. 
The accuracy of classification for the correct regime is denoted 
by the bold numbers, whereas the other percentages denote to 
what extent and where misclassifications occur. The accuracy 
of the classification schemes are evaluated using linear mea¬ 
surements (u in ( 1171) 1 with the cubic and quintic libraries illus¬ 
trated in Figs. 0 and [S] Also shown are classification results 
using nonlinear measurements (Q3 and U5 in 1181) . Nonlinear 
measurements, if possible, offer significant accuracy improve¬ 
ment and robustness to noise. 

Table HVl Note the clear improvement (100% accuracy) 
in using nonlinear measurements for classification tasks. 
Thus if the noise is driven by the sensor itself, then non¬ 
linear measurements may be quite advantageous. 

VII. RECONSTRUCTION AND THE 

GALERKIN-POD APPROXIMATION 

The classification step of the last section identifies the 
dynamical regime of the complex system by using spar¬ 
sity promoting £1 optimization on the learned libraries. 
Once the correct /3j regime is determined, reconstruc¬ 
tion of the solution and a future state prediction can 
be achieved through the POD-Galerkin approximation. 
Specifically, once the dynamical regime 0j has been iden¬ 
tified, then a subset of modes l —» form the 

correct modal basis for a POD-Galerkin approximation. 

To be more precise, recall that only a limited number of 
measurements are made as in (1131) . But now u = 
where the vector c is now the projection onto the smaller 
set of library modes associated with a single 0j. Thus 
instead of ED. we now we have 

u = P$ L;ft c. (19) 

Unlike the classification step, we can now determine c 
by simply solving the above equation using a standard 


Moore-Penrose pseudo-inverse operator f [ 35 } so that 
c = (P<I>£ a.)fu, i.e. it solves for c by minimizing the 
£2 norm. With c determined, the reconstruction of the 
solution thus follows: 

u = $ iA (P$ LA .) t u (20) 

This is the reconstruction of the system given the sparse 
measurement vector u and a classification (3j. The POD- 
Galerkin approximation for the future state can then be 
accomplished by using © and with the DEIM algorithm 
for evaluating the nonlinearities 0. The initial condi¬ 
tion for the POD-Galerkin is given from 03). Thus as 
advocated in previous work [7), JJ, accurate classification 
is accomplished with £1 optimization (decoding) while 
the more standard £2 norm is used for reconstruction 
and POD-Galerkin projection (encoding). Figure [2] illus¬ 
trates the execution state outlined here for classification, 
reconstruction and projection. 

VIII. CONCLUSIONS AND OUTLOOK 

In conclusion, we advocate a general theoretical frame¬ 
work for complex systems whereby low-rank libraries rep¬ 
resenting the optimal modal basis are constructed, or 
learned, from snapshot sampling of the dynamics. In or¬ 
der to make model reduction methods such as POD com¬ 
putationally efficient, especially in evaluating the nonlin¬ 
ear terms of the governing equations, nonlinear libraries 
are also constructed during the learning stage. This al¬ 
lows for the application of the discrete empirical interpo¬ 
lation method which identifies a limited number of spatial 
sampling locations that can allow for reconstruction of 
the nonlinear terms in a low-dimensional manner. Such 
sparse sampling of the nonlinearity is directly related to 
compressive sensing strategies whereby a small number 
of sensors can be used to characterize the dynamics of the 
complex system. Indeed, the POD method, when com¬ 
bined with DEIM and compressive sensing, can (i) cor¬ 
rectly identifying the dynamical parameter regime, (ii) 
reconstruct the full state dynamics and (iii) produce a 
low-rank prediction of the future state of the complex 
system. All of these tasks are accomplished in a low¬ 
dimensional way, unlike standard POD-Galerkin models 
whose nonlinearities can prove to be computationally in¬ 
efficient. 

To be more precise about our learning algorithm for the 
complex system, We construct the library modes repre¬ 
senting the dynamics by the £ 2 -optimal proper orthogo¬ 
nal decomposition. Several libraries are constructed: one 
for linear snapshot measurements, one for each nonlinear 
term, and one which combines all the nonlinear terms 
together with their prescribed weightings. The DEIM 
algorithm then allows us to identify sparse measurement 
locations capable of both classifying the dynamics regime 
of the complex system and efficiently evaluating the non¬ 
linear inner products for a POD-Galerkin projection of 
the system. Indeed, the dynamical state is identified from 
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limited noisy measurements using the sparsity promoting 
l\ norm and the compressive sensing architecture. The 
strategy for building modal libraries by concatenating 
truncated POD libraries across a range of relevant bifur¬ 
cation parameters may be viewed as a simple machine 
learning implementation. The resulting modal libraries 
are a natural sparse basis for the application of com¬ 
pressive sensing. After the expensive one-time library- 
building procedure, accurate identification, projection, 
and reconstruction may be performed entirely in a low¬ 
dimensional framework. 

With three DEIM determined sensor locations, it is 
possible to accurately classify bifurcation regimes, re¬ 
construct the low-dimensional content, and simulate the 
Galerkin projected dynamics of the complex Ginzburg 
Landau equation. In addition, we investigate the per¬ 
formance of sparse representation with the addition of 
sensor noise. For moderate noise levels, the method 
accurately classifies the correct dynamic regime. Non¬ 
linear measurements dramatically improve the classifica¬ 
tion procedure. Interestingly, the DIEMs algorithm not 
only provides nearly optimal sensor positioning, it also 
helps perform POD-Galerking truncations in a fully low- 
rank manner, thus avoiding the computational expense of 
evaluating nonlinear terms using the POD methodology. 


Overall, the combination of 1 2 low-rank representations 
and £1 sparse sampling enables efficient characterization 
and manipulation of low-rank dynamical systems. 

For modern complex systems, it is known that non¬ 
linearity plays a dominant role and shapes the under¬ 
lying spatio-temporal dynamics and modal structures, 
thus necessitating a new approach, such as that pre¬ 
sented here, for extracting these critical structures. As 
has been demonstrated, although nonlinearity drives new 
modal structures, it does not destroy the underlying low¬ 
dimensional nature of the dynamics. Methods that take 
advantage of such underlying structure are critical for 
developing theoretical understanding and garnering in¬ 
sight into the fundamental interactions of a vast array of 
physical, engineering and biological systems. 
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